This code generates the useful results for the article Guisset, Martin, and Govaerts (2019) (Cf. the Reference section below).

1 Load libraries and functions

Different librairies are needed to run this code.

Specify where are located the functions and data

# Choix du répertoire de travail

# directoryMain <- file.path('/Users/manon/Desktop/CodeArticleGuissetetAl')
# directoryMain <-
# file.path('~/Dropbox/PartageDiversProfessionel/Partage_Package_ASCA/CodeGuissetetAlMinimum')
directoryMain <- file.path("C:/Users/franc/Desktop/Stage_Package/RPackage")
directoryData <- file.path(paste0(directoryMain, "/Data"))

# lecture des routines R
DirectoryFun <- file.path(paste0(directoryMain, "/LIB"))
source(file.path(DirectoryFun, "Decomposition_functions.R"))
# source(file.path(DirectoryFun,'AComDim_functions.R'))
# source(file.path(DirectoryFun,'PARAFAC_array.R'))
# source(file.path(DirectoryFun,'AMOPLS_functions.R'))
source(file.path(DirectoryFun, "permutationTest.R"))
source(file.path(DirectoryFun, "matrixDecomposition.R"))


# output path
directoryOut <- file.path(paste0(directoryMain, "/Outputs"))

nperm <- 1000

2 Data

2.1 Uploading data

D_UCH.RData contains the complete design of the urine database and X_UCH.RData the complete spectral matrix. matrixDecomposition.R decomposes the spectral matrix according to the design

load(file.path(directoryData, "AllData.Rdata"))
pander("outcomes")

outcomes

str(head(outcomes))
##  num [1:6, 1:600] 0.0312 0.0581 0.027 0.0341 0.0406 ...
##  - attr(*, "dimnames")=List of 2
##   ..$   : chr [1:6] "M2C00D2R1" "M2C00D2R2" "M2C02D2R1" "M2C02D2R2" ...
##   ..$ X1: chr [1:600] "9.9917004" "9.9753204" "9.9590624" "9.9427436" ...
pander("design")

design

str(design)
## 'data.frame':    34 obs. of  5 variables:
##  $ Hippurate: Factor w/ 3 levels "0","1","2": 1 1 1 1 1 1 2 2 2 2 ...
##  $ Citrate  : Factor w/ 3 levels "0","2","4": 1 1 2 2 3 3 1 1 2 2 ...
##  $ Dilution : Factor w/ 1 level "diluted": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Day      : Factor w/ 2 levels "2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Time     : Factor w/ 2 levels "1","2": 1 2 1 2 1 2 1 2 1 2 ...
pander("formula")

formula

print(formula)
## outcomes ~ Hippurate + Citrate + Time + Hippurate * Citrate + 
##     Time * Hippurate + Time * Citrate + Hippurate * Citrate * 
##     Time
attach(design)

n <- dim(outcomes)[1]
m <- dim(outcomes)[2]

pander("table(Hippurate, Citrate, Time)")

table(Hippurate, Citrate, Time)

table(Hippurate, Citrate, Time)
## , , Time = 1
## 
##          Citrate
## Hippurate 0 2 4
##         0 2 1 1
##         1 2 2 2
##         2 2 2 2
## 
## , , Time = 2
## 
##          Citrate
## Hippurate 0 2 4
##         0 2 2 2
##         1 2 2 2
##         2 2 2 2
Variables = as.numeric(dimnames(outcomes)[[2]])  # Vector with variable names 
Samples = dimnames(outcomes)[[1]]  # Vector with sample names

2.2 Figure 1 spectrum

### Spectral profile
plot(outcomes["M2C24D2R1", ], type = "l", xaxt = "n", ylab = "Intensity", xlab = "ppm", 
    cex.lab = 1, mgp = c(2.3, 1, 0))
axis(side = 1, at = seq(1, m, 55), labels = round(as.numeric(colnames(outcomes))[seq(1, 
    m, 55)], 1))
title(expression(paste("Example of ", a^{
    1
}, "H NMR spectral profile")), line = 1)

3 Define a nice Scatterplot Matrix function

#------- define the Scatterplot Matrix function
# if spotPoints == TRUE, will circle specific outlying observations
pairsBG = function(MatScores, titre = "Scores", spotPoints = FALSE, ...) {
    panelup = function(x, y) {
        col <- rep("blue", length(Citrate))
        col[Citrate == "2"] <- "forestgreen"
        col[Citrate == "4"] <- "red"
        pch <- rep(4, length(Hippurate))
        pch[Hippurate == "1"] <- 16
        pch[Hippurate == "2"] <- 2
        
        points(x, y, col = col, pch = pch)
        
        if (spotPoints) {
            text(x, y, textlabel, pos = c(3), col = "darkturquoise")
            points(x, y, pch = pch2, col = "darkturquoise")
        }
        
    }
    paneldown = function(x, y) {
        
        pch = rep(1, length(Time))
        pch[Time == 2] <- 3
        
        col <- rep("orange", length(Time))
        col[Time == 2] <- "black"
        
        points(x, y, col = col, pch = pch)
    }
    pairs(MatScores, main = titre, upper.panel = panelup, lower.panel = paneldown, 
        gap = 0.3, ...)
}


#------- define a function to draw the legends next to the scatterplot

legendsScatterMatrix <- function() {
    legend("topright", title = expression(bold("Hippurate")), legend = c("0", 
        "1", "2"), bty = "n", col = c(1, 1), pch = c(4, 16, 2), inset = c(0.11, 
        0.1), cex = 0.7)
    
    legend("topright", title = expression(bold("Citrate")), legend = c("0", 
        "2", "4"), bty = "n", seg.len = 0.4, pch = 15, col = c("blue", "forestgreen", 
        "red"), inset = c(0.11, 0.4), cex = 0.7)
    
    legend("topright", title = expression(bold("Time")), legend = c("1", "2"), 
        bty = "n", seg.len = 0.4, pch = c(1, 3), col = c("orange", "black"), 
        inset = c(0.11, 0.7), cex = 0.7)
}

4 PCA on the spectral matrix

PCA decomposition is first applied on the spectral matrix outcomes by applying the function SVDforPCA from the package MBXUCL. Scores plots and loadings plots are obtained with DrawScores and DrawLoadings functions respectively.

4.1 PCA and % var explained

pcaOutcomes = SVDforPCA(outcomes)
eig.res = rbind(pcaOutcomes$var, pcaOutcomes$var * 100/sum(pcaOutcomes$var), 
    pcaOutcomes$cumvar)[, 1:6]
rownames(eig.res) = c("Variances", "Prop Var", "Cum Eigen Values")
pander(eig.res)
  PC1 PC2 PC3 PC4 PC5 PC6
Variances 45.25 24.83 16.75 6.716 3.412 0.9151
Prop Var 45.25 24.83 16.75 6.716 3.412 0.9151
Cum Eigen Values 45.25 70.08 86.82 93.54 96.95 97.87

4.2 scores plot

col <- c(`0` = "blue", `2` = "forestgreen", `4` = "red")
pch <- c(`0` = 4, `1` = 1, `2` = 2)

responsematrix_Scores <- DrawScores(obj = pcaOutcomes, type.obj = "PCA", drawNames = F, 
    createWindow = FALSE, main = "Reponse matrix scores plot", color = Citrate, 
    pch = Hippurate, size = 2, cex.lab = 3, axes = c(1, 2), xlab = NULL, ylab = NULL, 
    drawEllipses = FALSE, typeEl = "norm", levelEl = 0.9, drawPolygon = FALSE, 
    noLegend = FALSE, legend_color_manual = col, legend_shape_manual = pch) + 
    theme(panel.background = element_rect(fill = "white", colour = "grey50")) + 
    geom_text(aes(label = Time), hjust = c(0), vjust = -0.9, col = "black", 
        cex = 3)

responsematrix_Scores

Time <- paste0("time ", Time)

DrawScores(obj = pcaOutcomes, type.obj = "PCA", drawNames = F, createWindow = FALSE, 
    main = "Reponse matrix scores plot", color = Citrate, pch = Time, size = 2, 
    cex.lab = 3, axes = c(1, 2), xlab = NULL, ylab = NULL, drawEllipses = FALSE, 
    typeEl = "norm", levelEl = 0.9, drawPolygon = FALSE, noLegend = FALSE, legend_color_manual = col) + 
    scale_shape_manual(name = "Time", breaks = c("time 1", "time 2"), values = c("1", 
        "2"))

Time <- design$Time

DrawScores(pcaOutcomes, type.obj = "PCA", drawNames = TRUE, createWindow = F, 
    main = "Reponse matrix score plot", color = Citrate, pch = Hippurate, axes = c(1, 
        2), size = 2.5)

DrawScores(pcaOutcomes, type.obj = "PCA", drawNames = TRUE, createWindow = F, 
    main = "Reponse matrix score plot", color = Citrate, pch = Hippurate, axes = c(3, 
        4), size = 2.5)

DrawScores(pcaOutcomes, type.obj = "PCA", drawNames = TRUE, createWindow = F, 
    main = "Reponse matrix score plot", color = Citrate, pch = Hippurate, axes = c(5, 
        6), size = 2.5)

4.3 scatter plot des scores

par(xpd = TRUE)
pairsBG(pcaOutcomes$scores[, 1:6], "Scores PCA")
legendsScatterMatrix()

5 GLM decomposition

5.1 Model definition

ModelTerms_sansE <- attr(terms(formula), "term.labels")  # without residuals
ModelTerms <- c(ModelTerms_sansE, "residuals")  # with residuals

ModelTerms_abbrev <- ModelTerms  # Abbreviated model terms
index <- gregexpr(":", ModelTerms)
for (i in 1:length(ModelTerms)) {
    if (index[[i]][1] != -1) {
        ModelTerms_abbrev[i] <- substr(ModelTerms[i], 1, 1)
        index2 <- index[[i]]
        for (k in 1:length(index2)) {
            ModelTerms_abbrev[i] <- paste(ModelTerms_abbrev[i], substr(ModelTerms[i], 
                index2[k] + 1, index2[k] + 1), sep = "x")
        }
    }
}

ModelTerms_abbrev[length(ModelTerms_abbrev)] = "Residuals"

p <- length(ModelTerms)  # N model terms

5.2 GLM decomposition

# Décomposition de Y en matrice des effets par GLM

resGLM <- matrixDecomposition(formula, outcomes, design)
modelMatrix <- sapply(resGLM$modelMatrixByEffect, function(x) x)
nparam <- sum(sapply(modelMatrix, function(x) dim(x)[2])) - 1


# construction de la liste des matrices des effets purs
EffectMatGLM <- resGLM$effectMatrices[-1]  # minus intercept
res <- vector(mode = "list")
res[[1]] <- resGLM$residuals
EffectMatGLM <- c(EffectMatGLM, residuals = res)  # plus residuals
pander("names(EffectMatGLM)")

names(EffectMatGLM)

names(EffectMatGLM)
## [1] "Hippurate"              "Citrate"               
## [3] "Time"                   "Hippurate:Citrate"     
## [5] "Hippurate:Time"         "Citrate:Time"          
## [7] "Hippurate:Citrate:Time" "residuals"
# Calcul des matrices augmentées
EffectMatGLMAug <- resGLM$effectMatrices[-1]  # effectMatrices minus intercept
EffectMatGLMAug <- lapply(EffectMatGLMAug, function(x) x + resGLM$residuals)
EffectMatGLMAug <- c(EffectMatGLMAug, residuals = res)  # plus residuals
pander("names(EffectMatGLMAug)")

names(EffectMatGLMAug)

names(EffectMatGLMAug)
## [1] "Hippurate"              "Citrate"               
## [3] "Time"                   "Hippurate:Citrate"     
## [5] "Hippurate:Time"         "Citrate:Time"          
## [7] "Hippurate:Citrate:Time" "residuals"
Xmat <- resGLM$modelMatrix[, -1]  # modelMatrix minus intercept

5.3 Variation percentages

How to get effects matrices for all the effects. The sum is equal to 100% for balanced designs

VariationPercentages <- unlist(resGLM$variationPercentages)
pander(round(VariationPercentages, 2))
Table continues below
Hippurate Citrate Time Hippurate:Citrate Hippurate:Time
39.31 29.91 16.24 1.54 6.23
Citrate:Time Hippurate:Citrate:Time residuals
0.54 1.68 4.3 
barplot(t(unlist(resGLM$variationPercentages)), las = 2)

sum(VariationPercentages)
## [1] 99.74193

5.4 Permutation tests

5.4.1 Result permut

load(file.path(directoryOut, "PermutationResASCA.RData"))

names(pval) <- names(ASCAFSTAT) <- names(ASCAFSTATperm) <- ModelTerms_sansE

pander(pval)
Table continues below
Hippurate Citrate Time Hippurate:Citrate Hippurate:Time Citrate:Time
0 0 0 0.146 0 0.448
Hippurate:Citrate:Time
0.104 
for (i in ModelTerms_sansE) {
    hist(ASCAFSTATperm[[i]], xlim = c(0, max(ASCAFSTAT[i], ASCAFSTATperm[[i]])), 
        breaks = 20, main = i, freq = FALSE)
    points(ASCAFSTAT[i], 0, col = "red", pch = 19)
    rangex = range(ASCAFSTAT[i], ASCAFSTATperm[[i]])
    xf = seq(0, rangex[2], by = (rangex[2] - rangex[1])/1000)
    yf = df(xf, n - 1, n - 1)
    lines(xf, yf)
}

6 ASCA-GLM

6.1 scores/loadings

listgraphs <- list()
varASCA <- list()

for (i in 1:length(ModelTerms)) {
    # PCA on the non augmented effect matrices
    ascaSVD = SVDforPCA(EffectMatGLM[[i]])
    ascaSVD$scores = round(ascaSVD$scores, 5)
    varASCA[[i]] <- ascaSVD$var
    # Scores plots with Citrate/hippurate colors
    listgraphs[[paste0(ModelTerms[i], "-CH")]] <- DrawScores(ascaSVD, type.obj = "PCA", 
        drawNames = TRUE, createWindow = F, main = paste0(ModelTerms[i], "-CH", 
            " 
                                                                     scores plot - ASCA-GLM"), 
        color = as.factor(Citrate), pch = as.factor(Hippurate), axes = c(1, 
            2), size = 2.5)
    # Scores plots with Time/Day colors
    listgraphs[[paste0(ModelTerms[i], "-DM")]] <- DrawScores(ascaSVD, type.obj = "PCA", 
        drawNames = TRUE, createWindow = F, main = paste0(ModelTerms[i], "-DM", 
            " scores plot - ASCA-GLM"), color = as.factor(Day), pch = as.factor(Time), 
        axes = c(1, 2), size = 2.5)
    
    # Loadings plots
    listgraphs[[paste0(ModelTerms[i], "-Loadings")]] <- DrawLoadings(ascaSVD, 
        type.obj = "PCA", createWindow = F, main = paste0(ModelTerms[i], " loadings plot - ASCA-GLM"), 
        axes = c(1:2), loadingstype = "s", num.stacked = 2, xlab = NULL, ylab = NULL, 
        ang = "0", xaxis_type = "numerical", nxaxis = 10)
}

listgraphs
## $`Hippurate-CH`

## 
## $`Hippurate-DM`

## 
## $`Hippurate-Loadings`
## $`Hippurate-Loadings`[[1]]

## 
## 
## $`Citrate-CH`

## 
## $`Citrate-DM`

## 
## $`Citrate-Loadings`
## $`Citrate-Loadings`[[1]]

## 
## 
## $`Time-CH`

## 
## $`Time-DM`

## 
## $`Time-Loadings`
## $`Time-Loadings`[[1]]

## 
## 
## $`Hippurate:Citrate-CH`

## 
## $`Hippurate:Citrate-DM`

## 
## $`Hippurate:Citrate-Loadings`
## $`Hippurate:Citrate-Loadings`[[1]]

## 
## 
## $`Hippurate:Time-CH`

## 
## $`Hippurate:Time-DM`

## 
## $`Hippurate:Time-Loadings`
## $`Hippurate:Time-Loadings`[[1]]

## 
## 
## $`Citrate:Time-CH`

## 
## $`Citrate:Time-DM`

## 
## $`Citrate:Time-Loadings`
## $`Citrate:Time-Loadings`[[1]]

## 
## 
## $`Hippurate:Citrate:Time-CH`

## 
## $`Hippurate:Citrate:Time-DM`

## 
## $`Hippurate:Citrate:Time-Loadings`
## $`Hippurate:Citrate:Time-Loadings`[[1]]

## 
## 
## $`residuals-CH`

## 
## $`residuals-DM`

## 
## $`residuals-Loadings`
## $`residuals-Loadings`[[1]]

6.2 Variance decomposition for PC1 and PC2

names(varASCA) <- ModelTerms

pander("ASCA variance decomposition for PC1 and 2")

ASCA variance decomposition for PC1 and 2

# PC1
round(sapply(varASCA, function(x) x[1]), 2)
##              Hippurate.PC1                Citrate.PC1 
##                      97.71                      98.22 
##                   Time.PC1      Hippurate:Citrate.PC1 
##                     100.00                      44.01 
##         Hippurate:Time.PC1           Citrate:Time.PC1 
##                      93.92                      90.76 
## Hippurate:Citrate:Time.PC1              residuals.PC1 
##                      47.23                      48.54
# PC2
round(sapply(varASCA, function(x) x[2]), 2)
##              Hippurate.PC2                Citrate.PC2 
##                       2.29                       1.78 
##                   Time.PC2      Hippurate:Citrate.PC2 
##                       0.00                      38.51 
##         Hippurate:Time.PC2           Citrate:Time.PC2 
##                       6.08                       9.24 
## Hippurate:Citrate:Time.PC2              residuals.PC2 
##                      27.49                      16.90
varASCAPC12 <- rbind(sapply(varASCA, function(x) x[1]), sapply(varASCA, function(x) x[2]))


pander("sum over PC1 and PC2")

sum over PC1 and PC2

round(apply(varASCAPC12, 2, sum), 2)
##              Hippurate.PC1                Citrate.PC1 
##                     100.00                     100.00 
##                   Time.PC1      Hippurate:Citrate.PC1 
##                     100.00                      82.52 
##         Hippurate:Time.PC1           Citrate:Time.PC1 
##                     100.00                     100.00 
## Hippurate:Citrate:Time.PC1              residuals.PC1 
##                      74.72                      65.44

6.3 Hippurate scores and loadings

# Hippurate scores and loadings
ascaSVD <- SVDforPCA(EffectMatGLM$Hippurate)
ascaSVD$scores <- round(ascaSVD$scores, 5)

ASCAScoresHippurate <- DrawScores(ascaSVD, type.obj = "PCA", drawNames = F, 
    createWindow = F, main = "ASCA scores plot - Hippurate effect", axes = c(1, 
        2), size = 2.5, pch = Hippurate, color = Hippurate, noLegend = TRUE) + 
    scale_colour_manual(name = "Hippurate:", breaks = c("0", "1", "2"), values = c("blue", 
        "forestgreen", "red")) + scale_shape_manual(name = "Hippurate:", breaks = c("0", 
    "1", "2"), values = c(4, 1, 2))


ASCAScoresHippurate

ASCALoadingHippurate <- DrawLoadings(ascaSVD, type.obj = "PCA", createWindow = F, 
    main = "ASCA loading plot - Hippurate effect", axes = c(1), loadingstype = "s", 
    num.stacked = 2, xlab = "ppm", ylab = "", ang = "0", xaxis_type = "numerical", 
    nxaxis = 10)

ASCALoadingHippurate
## [[1]]

7 ASCA-E-GLM

7.1 Overview of scores and loadings

PCA decomposition is applied on each non-augmented spectral matrix. Residuals are then projected on the dimensions 2 first dimension for each effect et scores plots prepared. The loadings are just those of the original PCA.

listgraphs <- list()
for (i in 1:(length(ModelTerms))) {
    # PCA on the non-augmented effect matrices
    ascaSVD <- SVDforPCA(EffectMatGLM[[i]])
    # scores 1 and 2 correction by adding a projection of the residuals matrix
    ascaSVD$scores[, 1:2] <- (EffectMatGLM[[i]] + EffectMatGLM[[p]]) %*% ascaSVD$loadings[, 
        1:2]
    # Scores plot with Citrate/hippurate colors
    listgraphs[[paste0(ModelTerms[i], "-CH")]] <- DrawScores(ascaSVD, type.obj = "PCA", 
        drawNames = F, createWindow = F, main = paste0(ModelTerms[i], "-CH", 
            " score plot"), color = Citrate, pch = Hippurate, axes = c(1, 2), 
        size = 2.5)
    # Scores plot with Media/Time colors
    listgraphs[[paste0(ModelTerms[i], "-RM")]] <- DrawScores(ascaSVD, type.obj = "PCA", 
        drawNames = F, createWindow = F, main = paste0(ModelTerms[i], "-DM", 
            " score plot"), color = Time, pch = Dilution, axes = c(1, 2), size = 2.5)
    # Loadings plot
    listgraphs[[paste0(ModelTerms[i], "-Loadings")]] <- DrawLoadings(ascaSVD, 
        type.obj = "PCA", createWindow = F, main = paste0(ModelTerms[i], " loading plot"), 
        axes = c(1:2), loadingstype = "s", num.stacked = 2, xlab = NULL, ylab = NULL, 
        ang = "0", xaxis_type = "numerical", nxaxis = 10)
}
listgraphs
## $`Hippurate-CH`

## 
## $`Hippurate-RM`

## 
## $`Hippurate-Loadings`
## $`Hippurate-Loadings`[[1]]

## 
## 
## $`Citrate-CH`

## 
## $`Citrate-RM`

## 
## $`Citrate-Loadings`
## $`Citrate-Loadings`[[1]]

## 
## 
## $`Time-CH`

## 
## $`Time-RM`

## 
## $`Time-Loadings`
## $`Time-Loadings`[[1]]

## 
## 
## $`Hippurate:Citrate-CH`

## 
## $`Hippurate:Citrate-RM`

## 
## $`Hippurate:Citrate-Loadings`
## $`Hippurate:Citrate-Loadings`[[1]]

## 
## 
## $`Hippurate:Time-CH`

## 
## $`Hippurate:Time-RM`

## 
## $`Hippurate:Time-Loadings`
## $`Hippurate:Time-Loadings`[[1]]

## 
## 
## $`Citrate:Time-CH`

## 
## $`Citrate:Time-RM`

## 
## $`Citrate:Time-Loadings`
## $`Citrate:Time-Loadings`[[1]]

## 
## 
## $`Hippurate:Citrate:Time-CH`

## 
## $`Hippurate:Citrate:Time-RM`

## 
## $`Hippurate:Citrate:Time-Loadings`
## $`Hippurate:Citrate:Time-Loadings`[[1]]

## 
## 
## $`residuals-CH`

## 
## $`residuals-RM`

## 
## $`residuals-Loadings`
## $`residuals-Loadings`[[1]]

7.2 PCA on the effect matrices

# PC var for each effect matrix PCA
pcVar <- vector("list", length = length(ModelTerms))
Y_Loadings <- vector("list", length = length(ModelTerms))
ascaSVD_scores <- vector("list", length = (length(ModelTerms) + 1))
names(pcVar) <- names(Y_Loadings) <- ModelTerms

names(ascaSVD_scores) <- c(ModelTerms[-p], "residuals1", "residuals2")

for (i in 1:(p - 1)) {
    ascaSVD = SVDforPCA(EffectMatGLM[[i]])
    pcVar[[ModelTerms[i]]] <- ascaSVD$var/100
    Y_Loadings[[ModelTerms[i]]] <- ascaSVD$loadings
    ascaSVD$scores[, 1:2] = (EffectMatGLM[[i]] + EffectMatGLM$residuals) %*% 
        ascaSVD$loadings[, 1:2]
    ascaSVD_scores[[ModelTerms[i]]] <- ascaSVD$scores[, 1]
}

# Error
ascaSVD = SVDforPCA(EffectMatGLM$residuals)
ascaSVD_error <- ascaSVD
pcVar[["residuals"]] <- ascaSVD$var/100
ascaSVD$scores = round(ascaSVD$scores, 5)
ascaSVD_scores[["residuals1"]] <- ascaSVD$scores[, 1]
ascaSVD_scores[["residuals2"]] <- ascaSVD$scores[, 2]

7.3 Saliences plot

saliences <- matrix(0, nrow = p, ncol = p + 1)
colnames(saliences) <- paste(c(ModelTerms[-p], "residuals1", "residuals2"), 
    "PC1")

colnames(saliences)[p + 1] <- "residuals2 PC2"

rownames(saliences) <- ModelTerms
for (i in 1:(p - 1)) {
    saliences[ModelTerms[i], colnames(saliences)[i]] <- pcVar[[i]][1]
}

saliences["residuals", "residuals1 PC1"] <- pcVar$residuals[1]
saliences["residuals", "residuals2 PC2"] <- pcVar$residuals[2]


# plot --------------------------------------------------------------------

ylim <- c(0, 1)
angle1 <- rep(c(45, 45, 135), length.out = p)
angle2 <- rep(c(45, 135, 135), length.out = p)
density1 <- seq(5, 35, length.out = p)
density2 <- seq(5, 35, length.out = p)
col <- rainbow((p))


par(xpd = TRUE, mar = c(5, 4, 4, 4))
x <- barplot(saliences, beside = TRUE, ylim = ylim, xaxt = "n", col = col, angle = angle1, 
    density = density1, main = "Effects explained by ASCA-E components")
barplot(saliences, add = TRUE, beside = TRUE, ylim = ylim, xaxt = "n", col = col, 
    angle = angle2, density = density2)

labs <- as.vector(colnames(saliences))
text(cex = 1, x = colMeans(x) - 0.25, y = -0.1, labs, xpd = TRUE, srt = 45)

7.4 ASCA variance decomposition

pcexp_pereffect_Resid <- (ascaSVD_error$var[1:2]/100) * VariationPercentages["residuals"]

pcexp_pereffect_Effects <- VariationPercentages * round(sapply(varASCA, function(x) x[1]), 
    2)/100

pcexp_pereffect_perPC <- c(pcexp_pereffect_Effects[-p], pcexp_pereffect_Resid)

pander("pcexp_pereffect")

pcexp_pereffect

pander(pcexp_pereffect_perPC)
Table continues below
Hippurate Citrate Time Hippurate:Citrate Hippurate:Time
38.41 29.37 16.24 0.6789 5.85
Citrate:Time Hippurate:Citrate:Time PC1 PC2
0.4889 0.7954 2.086 0.7263

7.5 Scores plots

ascaSVD <- SVDforPCA(EffectMatGLM$Hippurate)
ascaSVD$scores[, 1:2] <- (EffectMatGLM$Hippurate + EffectMatGLM$residuals) %*% 
    ascaSVD$loadings[, 1:2]

col <- c(`0` = "blue", `1` = "forestgreen", `2` = "red")

ASCAEScoresHippurate <- DrawScores(ascaSVD, type.obj = "PCA", drawNames = F, 
    createWindow = F, main = "ASCA-E scores plot - Hippurate effect", axes = c(1, 
        2), size = 2.5, noLegend = FALSE, pch = Hippurate, color = Hippurate) + 
    scale_colour_manual(name = "Hippurate:", breaks = c("0", "1", "2"), values = c("blue", 
        "forestgreen", "red")) + scale_shape_manual(name = "Hippurate:", breaks = c("0", 
    "1", "2"), values = c(4, 1, 2))

ASCAEScoresHippurate

# Citrate scores
ascaSVD <- SVDforPCA(EffectMatGLM$Citrate)
ascaSVD$scores[, 1:2] <- (EffectMatGLM$Citrate + EffectMatGLM$residuals) %*% 
    ascaSVD$loadings[, 1:2]

col <- c(`0` = "black", `2` = "orange", `4` = "red")
DrawScores(ascaSVD, type.obj = "PCA", drawNames = F, createWindow = F, main = "ASCA-E scores plot - Citrate effect", 
    axes = c(1, 2), size = 2.5, noLegend = TRUE, pch = Citrate, color = Citrate) + 
    scale_colour_manual(name = "Citrate:", breaks = c("0", "2", "4"), values = c("black", 
        "orange", "red")) + scale_shape_manual(name = "Citrate:", breaks = c("0", 
    "2", "4"), values = c(25, 15, 4))

7.5.1 Scores scatter plots

# scores plot
mat_ascaSVD <- do.call(cbind, ascaSVD_scores)
names(ascaSVD_scores)
## [1] "Hippurate"              "Citrate"               
## [3] "Time"                   "Hippurate:Citrate"     
## [5] "Hippurate:Time"         "Citrate:Time"          
## [7] "Hippurate:Citrate:Time" "residuals1"            
## [9] "residuals2"
ascaSVD_scores <- c(ascaSVD_scores[p:(p + 1)], ascaSVD_scores[1:(p - 1)])

ascaSVD_scores$Time <- ascaSVD_scores$Time * -1
ascaSVD_scores$`Hippurate:Time` <- ascaSVD_scores$`Hippurate:Time` * -1


labels <- paste0(c(paste(ModelTerms_abbrev, "\n PC 1"), "Residuals \n PC 2"), 
    "\n ", round(pcexp_pereffect_perPC, 2), "%")

# plot --------------------------------------------------------------------
par(xpd = TRUE)
pairsBG(mat_ascaSVD, titre = "ASCA-E Scores plots", oma = c(3, 3, 5, 15), labels = labels, 
    spotPoints = FALSE)

legendsScatterMatrix()

7.6 Loadings plots

# loadings plot
loadingsAPCAGLM <- do.call(cbind, sapply(Y_Loadings, function(x) x[, 1]))

colnames(loadingsAPCAGLM) <- paste(colnames(loadingsAPCAGLM), "PC1")
colnames(loadingsAPCAGLM)[colnames(loadingsAPCAGLM) == "Hippurate:Time PC1"] <- "HxT PC1"

loadingsAPCAGLM[, "Time PC1"] <- loadingsAPCAGLM[, "Time PC1"] * -1
loadingsAPCAGLM[, "HxT PC1"] <- loadingsAPCAGLM[, "HxT PC1"] * -1

# plot --------------------------------------------------------------------
LinePlot(t(loadingsAPCAGLM[, c(1, 2, 3, 5)]), main = "ASCA-GLM Loadings", type = "s", 
    num.stacked = 4, xaxis_type = "numerical", nxaxis = 10)
## [[1]]

8 APCA-GLM scores and loadings

listgraphs <- list()
for (i in 1:length(ModelTerms)) {
    # PCA on the augmented effect matrices
    ascaSVD <- SVDforPCA(EffectMatGLMAug[[i]])
    ascaSVD$scores <- round(ascaSVD$scores, 5)
    # Scores plot with Citrate/hippurate color
    listgraphs[[paste0(ModelTerms[i], "-CH")]] <- DrawScores(ascaSVD, type.obj = "PCA", 
        drawNames = F, createWindow = F, main = paste0(ModelTerms[i], "-CH", 
            " score plot"), color = Citrate, pch = Hippurate, axes = c(1, 2), 
        size = 2.5)
    # Scores plot with Media/Day color
    listgraphs[[paste0(ModelTerms[i], "-DM")]] <- DrawScores(ascaSVD, type.obj = "PCA", 
        drawNames = F, createWindow = F, main = paste0(ModelTerms[i], "-DM", 
            " score plot"), color = Day, pch = Dilution, axes = c(1, 2), size = 2.5)
    # Loadings plot
    listgraphs[[paste0(ModelTerms[i], "-Loadings")]] <- DrawLoadings(ascaSVD, 
        type.obj = "PCA", createWindow = F, main = paste0(ModelTerms[i], " loading plot"), 
        axes = c(1:2), loadingstype = "s", num.stacked = 2, xlab = NULL, ylab = NULL, 
        ang = "0", xaxis_type = "numerical", nxaxis = 10)
}
listgraphs
## $`Hippurate-CH`

## 
## $`Hippurate-DM`

## 
## $`Hippurate-Loadings`
## $`Hippurate-Loadings`[[1]]

## 
## 
## $`Citrate-CH`

## 
## $`Citrate-DM`

## 
## $`Citrate-Loadings`
## $`Citrate-Loadings`[[1]]

## 
## 
## $`Time-CH`

## 
## $`Time-DM`

## 
## $`Time-Loadings`
## $`Time-Loadings`[[1]]

## 
## 
## $`Hippurate:Citrate-CH`

## 
## $`Hippurate:Citrate-DM`

## 
## $`Hippurate:Citrate-Loadings`
## $`Hippurate:Citrate-Loadings`[[1]]

## 
## 
## $`Hippurate:Time-CH`

## 
## $`Hippurate:Time-DM`

## 
## $`Hippurate:Time-Loadings`
## $`Hippurate:Time-Loadings`[[1]]

## 
## 
## $`Citrate:Time-CH`

## 
## $`Citrate:Time-DM`

## 
## $`Citrate:Time-Loadings`
## $`Citrate:Time-Loadings`[[1]]

## 
## 
## $`Hippurate:Citrate:Time-CH`

## 
## $`Hippurate:Citrate:Time-DM`

## 
## $`Hippurate:Citrate:Time-Loadings`
## $`Hippurate:Citrate:Time-Loadings`[[1]]

## 
## 
## $`residuals-CH`

## 
## $`residuals-DM`

## 
## $`residuals-Loadings`
## $`residuals-Loadings`[[1]]

# Hippurate scores and loadings
ascaSVD <- SVDforPCA(EffectMatGLMAug$Hippurate)
ascaSVD$scores <- round(ascaSVD$scores, 5)

APCAScoresHippurate <- DrawScores(ascaSVD, type.obj = "PCA", drawNames = F, 
    createWindow = F, main = "APCA scores plot - Hippurate effect", axes = c(1, 
        2), noLegend = TRUE, size = 2.5, pch = Hippurate, color = Hippurate) + 
    scale_colour_manual(name = "Hippurate:", breaks = c("0", "1", "2"), values = c("blue", 
        "forestgreen", "red")) + scale_shape_manual(name = "Hippurate:", breaks = c("0", 
    "1", "2"), values = c(4, 1, 2)) + guides(colour = guide_legend(override.aes = list(shape = c(25, 
    15, 4), color = c("blue", "red", "red"))))

APCAScoresHippurate

APCALoadingHippurate <- DrawLoadings(ascaSVD, type.obj = "PCA", createWindow = F, 
    main = "APCA loading plot - Hippurate effect", axes = c(1), loadingstype = "s", 
    num.stacked = 2, xlab = "ppm", ylab = "", ang = "0", xaxis_type = "numerical", 
    nxaxis = 10)

9 Hippurate scores and loadings for ASCA(-E) and APCA

scoresASCAhip <- plot_grid(ASCAScoresHippurate, APCAScoresHippurate, ASCAEScoresHippurate, 
    align = "none", nrow = 1, rel_widths = c(0.3, 0.32, 0.39))

loadingsASCAhip <- plot_grid(ASCALoadingHippurate[[1]] + ylim(-0.2, 0.5), APCALoadingHippurate[[1]] + 
    ylim(-0.2, 0.5), align = "none", nrow = 1, rel_heights = c(1/2, 1/2))


gridExtra::grid.arrange(scoresASCAhip, loadingsASCAhip, heights = c(1.5, 1))

References

Guisset, S., M. Martin, and B. Govaerts. 2019. “Comparison of PARAFASCA, AComDim, and AMOPLS Approaches in the Multivariate GLM Modelling of Multi-Factorial Designs.” Chemometrics and Intelligent Laboratory Systems 184 (January): 44–63. doi:10.1016/j.chemolab.2018.11.006.